##9/19/2015
##Replication_Code_7
## This script reproduces Table 2 
##-- main regression results
rm(list=ls())
setwd("./../data/") # C:/Users/as9934/Dropbox/HansardProject/burstiness/bjps_replication/data/")

#get the pooled data, formed as described in text
load("pooled.rdata")

require(apsrtable)

## data contains following variables --- 

# member -- member id
# cabinet_member -- binary, in cabinet or not
# outlier -- in burstiness terms
# session_number -- which session is it? (numeric, counting from 1832)
# burst -- burstiness of the member
# prior_service -- did member previously have cabinet office?
# speeches -- how many speeches did the member make?
# speech_outlier-- was member an outlier in speech terms?

#mod 1 is the baseline, no interaction
mod1 <- glm(cabinet_member ~ outlier + session_number, data=pool.data, family="binomial")

#mod 2 uses speeches as an alternate measure (fit is worse)
mod2 <- glm(cabinet_member ~ speeches + session_number, data=pool.data, family = 'binomial') 

#mod 3 allows the interaction between outlying and session number
mod3 <- glm(cabinet_member ~outlier*session_number, data=pool.data, family='binomial')

#mod 4 is mod 3, but controls for prior service
mod4 <- glm(cabinet_member ~ outlier*session_number + prior_service , data=pool.data, family='binomial')

## cluster the SEs
#use Drew Dimmery's code: http://drewdimmery.com/robust-ses-in-r/
robust.se <- function(model, cluster){
  require(sandwich)
  require(lmtest)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- model$rank
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  rcse.se <- coeftest(model, rcse.cov)
  return(list(rcse.cov, rcse.se))
}


#function to do clustered ses.
clus <- function(mod=mod1){
  #Create the new variable with appropriate level names
  # (Dimmery suggests this)
  clustervar<-mapply(paste,"MP.",pool.data$member,sep="")
  #Save the coefficient test output to an element in the model object
  # -- get rid of mod5 <-> mod
  mod$coeftest<-robust.se(mod,clustervar)[[2]]
  mod$se <- robust.se(mod,clustervar)[[2]][,2]
  mod
}


#main results table
print(
  apsrtable(clus(mod1), clus(mod2), clus(mod3), clus(mod4), digits=4, stars='default')
)


